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Constraining Dark Matter Models from a Combined Analysis of Milky Way Satellites 

with the Fermi Large Area Telescope 
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Satellite galaxies of the Milky Way are among the most promising targets for dark matter searches 
in gamma rays. We present a search for dark matter consisting of weakly interacting massive 
particles, applying a joint likelihood analysis to 10 satellite galaxies with 24 months of data of the 
Fermi Large Area Telescope. No dark matter signal is detected. Including the uncertainty in the dark 
matter distribution, robust upper limits are placed on dark matter annihilation Cross sections. The 
95% confidence level upper limits range from about 10~ 26 cmV 1 at 5 GeV to about 5x 1CT 23 cm 3 s^ 1 
at 1 TeV, depending on the dark matter annihilation final state. For the first time, using gamma 
rays, we are able to rule out models with the most generic cross section (~3x 10“ 6 cm 3 s -1 for a 
purely s-wave cross section), without assuming additional boost factors. 


INTRODUCTION 

It is well-established that baryons contribute only 
about 20% of the mass density of matter in the Uni- 
verse [1] . The nature of the remaining 80% of matter, 
known as dark matter (DM), remains a mystery. One 
leading candidate consists of weakly interacting massive 


particles (WIMP), predicted in several extensions of the 
standard model of particle physics. If the WIMP is a Ma- 
jorana fermion, its pair annihilation will produce gamma 
rays with a flux given by = (cr aRn v) / (8nm^r) x 

Nw(E) x J(ip), where (cr atm v) is the velocity averaged 
pair annihilation cross section, m \ y is the WIMP mass, 
Nw(E) is the gamma-rav energy distribution per anni- 
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hilation. and J(ip) — j] An dldClp 2 [l(ip)\ is the line- 
of-sight (l.o.s.) integral of the squared DM density, p, 
toward a direction of observation, ip, integrated over a 
solid angle, A ft (see e.g., [2]; see also [3] for a review). 

Regions of local DM density enhancements with large 
■J(ip), or J factors, are potentially good targets for DM 
searches. Dwarf spheroidal satellite galaxies (dSphs) of 
the Milky Way are DM-dominated systems without ac- 
tive star formation or detected gas content [4, 5]. Thus, 
though the expected number of signal counts is not as 
high as from the Galactic center for instance, dSphs ex- 
hibit a favorable signal to noise ratio, and upper limits 
on a gamma-ray signal from DM annihilation have been 
obtained by the Fermi Large Area Telescope ( Fermi - 
LAT) [6, 7] as well as air Cherenkov telescopes [8-11]. 

In this Letter, we present new Fermi - LAT results on 
dSphs, with an updated dataset and two significant im- 
provements over our previous analyses: first, we com- 
bine all the dSph observations into a single joint likeli- 
hood function, which improves the statistical power of 
the analysis, and second, we take into account the uncer- 
tainties in estimates of the J-factors, thereby making our 
results more robust. 


FERMI-LAT OBSERVATIONS 

The Ferroi-LAT. the main instrument on board the 
Fermi observatory, is a pair-conversion telescope that 
detects gamma rays in the energy range from 20 MeV to 
> 300 GeV with unprecedented sensitivity. Further de- 
tails on the instrument can be found in [12], and current 
official performance figures are available in [13] . 

In this Letter, we use 24 months of Fermi - LAT data, 
recorded between 2008-08-04 and 2010-08-04, and the 
data reduction is performed with the Fermi-LAT data 
analysis package, ScienceTools [14]. Only “diffuse” class 
events with energy between 200 MeV and 100 GeV are 
used. To avoid contamination from Earth limb gamma 
rays, events with zenith angles larger than 100° are re- 
jected and time intervals when the observed sky position 
is occulted by the Earth are discarded from the lifetime 
calculation. We extract from this dataset regions of in- 
terest (ROIs) of radius 10° around the position of each 
dSph specified in Table I. 

In this Letter we add Segue 1 and Carina to the sample 
of 8 dSphs analyzed in [6], where further details on the 
selection criteria are provided. Carina has been added as 
two years of data now allow us to reasonably model the 
Galactic diffuse background at Galactic latitudes about 
-22°. Segue 1 has been added because there has been 
significant progress in estimating its DM distribution re- 
cently [18-20]; however, uncertainties in these estimates 
remain large. 

This Letter uses the instrument response functions 
P6V3 1 12] for the “diffuse” class of events. We also use the 


TABLE I. Position, distance, and J factor (under assumption 
of a Navarro-Frenk-White profile) of each dSph. The 4th col- 
umn shows the mode of the posterior distribution of log 10 J, 
and the 5th column indicates its 68% C.L. error. See the text 
for further details. The J factors correspond to the pair anni- 
hilation flux coming from a cone of solid angle AO = 2.4- ICrt 4 
sr. The final column indicates the reference for the kinematic 
dataset used. 


Name 

1 

deg. 

b 

deg. 

d log 10 (J) o- 

kpc log ln [GeV 2 cm~ 5 ] 

ref. 

Bootes I 

358.08 

69.62 

60 

17.7 

0.34 

[15] 

Carina 

260.11 

-22.22 

101 

18.0 

0.13 

[16] 

Coma Berenices 

241.9 

83.6 

44 

19.0 

0.37 

[17] 

Draco 

86.37 

34.72 

80 

18.8 

0.13 

[16] 

Fornax 

237.1 

-65.7 

138 

17.7 

0.23 

[16] 

Sculptor 

287.15 

-83.16 

80 

18.4 

0.13 

[16] 

Segue 1 

220.48 

50.42 

23 

19.6 

0.53 

[18] 

Sextans 

243.4 

42.2 

86 

17.8 

0.23 

[16] 

Ursa Major II 

152.46 

37.44 

32 

19.6 

0.40 

[17] 

Ursa Minor 

104.95 

44.80 

66 

18.5 

0.18 

[16] 


diffuse emission model derived and recommended by the 
Fermi - LAT Collaboration [21]. It includes the Galactic 
diffuse emission component ( glLiem-.v02.ftt ), and a cor- 
responding isotropic component ( isotropicjiem.jv02.txt ) 
that accounts for isotropic background light, unresolved 
sources and residual cosmic-ray contamination. Point 
sources from the 1FGL catalog [22] within 15° of each 
dSph (and a few additional faint sources detected in two 
years of data) are included in the model. A potential DM 
signal in each ROI is modeled as a point source wffiere the 
gamma-ray yields are obtained from the DMFit package 
[23] based on DarkSUSY [24], as implemented in the Sci- 
enceTools. For the J factors (defined in the Introduc- 
tion), we use the updated values summarized in Table I, 
which were estimated as described in the next section. 


J FACTORS FROM STELLAR VELOCITY DATA 

J factors are calculated using the line-of-sight velocities 
of the stars in the dSph and the Jeans equation via a 
Bayesian method as described in the literature (e.g., [25- 
27]). The mass of DM within the half-light radii of the 
dSphs is well-constrained, independent of the assumption 
of whether there is a core or a cusp in the central DM 
density distribution [16, 28]. 

We assume that the inner DM density profile scales 
as 1/r, a close match to the results seen in dark-matter- 
only simulations where the particles are initially cold (like 
WIMPs). Baryonic processes may alter the density pro- 
files in dSphs [29] . However, present velocity data are un- 
able to differentiate between cores and cusps in dSphs in 
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a model-independent manner. If the dSphs have constant 
density cores, then, in order to match the stellar velocity 
data constraint (essentially the dynamical mass within 
the half-light radius), the normalization of the density 
profile at the half-light radius would have to be increased 
(compared to the 1/r profile). For large constant density 
cores (comparable to or larger than the dSph half-light 
radius), this results in a larger J factor if the pair anni- 
hilation flux is integrated over a solid angle larger than 
that encompassing half the stellar luminosity. This is due 
to the fact that flux is dominated by annihilations in the 
outer parts for 1/r and shallower dark matter density 
profiles. For small cores, the J factor can be smaller but 
the change is proportionally smaller also. 

The observed half-light radii of the dSphs is less than 
or close to 0.5° (which is the radius corresponding to the 
solid angle of Afi = 2.4 • 1CP 4 sr used to compute the J 
factor). Thus, if we were to adopt a cored dark matter 
profile, the J factors for most of the dSphs would either 
increase or not change much. We have not attempted to 
model the possible correlations in the J-factor estimates 
of the different dSphs that would arise due to common 
baryonic feedback processes in these systems. These pro- 
cesses could, for example, create large constant density 
cores in the dark matter halos of all the dSphs. With a 
deeper understanding of galaxy formation on these small 
scales, it may be possible to refine the present constraint. 

The DM mass distribution as a function of the 
radius from the center of the dwarf is modeled 
as a Navarro-Frenk- White profile given by p(r) — 
0.08 V"max r s /[Gr(r + r,) 2 ], where V max is the maximum 
circular velocity possible for the dark matter halo. 
For this profile, the J factor in units of GeV 2 cm~ 5 
~ 1 0 l ' ( V max / 1 0 krri s“ 1 ) 4 (k pc/r, s ) ( 100 kpc/r/) 2 up to a 
function of (d/r s )(AQ./ir) 1 / 2 that is of the order of unity 
for parameters of interest. 

The stellar velocities used in the calculations are taken 
from the references listed in Table I. For the 6 classical 
dwarfs, we used the available velocity dispersion data in 
radial bins [16], and, for the fainter dwarfs (discovered in 
SDSS), we used the individual stellar velocities [15, 171 . 
We used a Gaussian distribution for the line-of-sight ve- 
locity measurements, adding intrinsic velocity dispersion 
and measurement error in quadrature (see Eq. 13 of [27]) 
and imposed spherical symmetry. For the binned velocity 
dispersion data, we used an approximation starting with 
the same Gaussian distribution for velocities and then 
assuming that the intrinsic velocity dispersion dominates 
the average measurement error (see Eq. 14 of [27]). From 
tests on a few dwarfs, we expect that this approximation 
could introduce a bias of about 50% to the most probable 
individual J factors. Other approximate likelihoods we 
tested also resulted in similar biases compared to Gaus- 
sian distribution for velocities. We assume a flat prior 
in ln(V max ) and a prior for r s given V max consistent with 
both Aquarius and Via Lactea II simulations [30, 31]. 


For the dSphs with the highest-quality datasets (i.e. the 
ones with the most stars and smallest errors, including 
Draco and Ursa Minor), the results do not change sig- 
nificantly if the flat ln(V' raax ) prior is changed to match 
the V m ax distribution of subhalos in A cold dark mat- 
ter simulations [27] or if we assume a flat ln(r s ) prior. 
However, for dSphs with sparse datasets, such as ultra- 
faint Ursa Major II (20 stars) and Segue 1 (66 stars), the 
results are prior dependent. For example, adopting the 
subhalo prior for V max decreases the median J by a fac- 
tor of - 2 for Segue 1 and Ursa Major II. The ultrafaint 
dSphs are promising candidates, but these and other sig- 
nificant uncertainties remain in the estimates of their DM 
halo mass. Considerable progress in dealing with some 
of these uncertainties has been made for Segue 1 [18-20], 
but we have opted to treat both Segue 1 and Ursa Major 
II in the same fashion as the other dSphs for the sake of 
uniformity in treating the priors. This is a limitation of 
the analysis at present, so we quote constraints with and 
without Segue 1 and Ursa Major II below. The final re- 
sults for the J factors within Afl = 2.4- 10~ 4 sr are listed 
in Table I. To be conservative, we assume no contribu- 
tion to the flux from DM substructure in the dSphs. The 
posterior distribution as well as the likelihood function 
for J are -well-described by a log normal function, which 
is used in order to include the uncertainty on J in the 
confidence interval calculation, as described in the next 
section. 


DATA ANALYSIS 


The ScienceTools analysis package is used to perform 
a binned Poisson likelihood fit to both spatial and spec- 
tral information in the data, with 30 energy bins loga- 
rithmically spaced from 200 MeV to 100 GeV and 10° 
square spatial maps with a bin size of 0.1°. The nor- 
malizations of the two diffuse components are left free in 
all EOIs, together with the normalizations of the point 
sources within 5° of the dSph position. The first improve- 
ment to the analysis in [6] consists of combining the DM 
signal across all the ROIs. Indeed, the J factor is differ- 
ent for each dSph, but the characteristics of the WIMP 
candidate {raw-, {c ant iu), annihilation channels and their 
branching ratios) can be assumed to be universal. As 
a consequence, the Ferrm-LAT Collaboration developed 
the CompositeS code in the ScienceTools, to allow tying 
any set of parameters across any set of ROIs. The sec- 
ond improvement is that uncertainties on the J factor are 
taken into account in the fit procedure by adding another 
term to the likelihood that represents the measurement 
uncertainties. With this addition, the joint, likelihood 
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considered in our analysis becomes 


L(D|pw,{p}i) = n^PIPw.Pi 


In (10) Jj\/2TTaj 


- [*Og 1 0 ( A ) - lo S 10 ( A ) j 2 / 2cr i 


CD 


where LJ jAT denotes the binned Poisson likelihood that is 
commonly used in a standard single ROI analysis of the 
L AT data and takes full account of the point-spread func- 
tion, including its energy dependence; i indexes the ROIs; 
D represents the binned gamma-ray data; pw represents 
the set of ROI-independent DM parameters ((<y ann f) and 
rriw)\ and {p} t are the ROI-dependent model parame- 
ters. In this analysis, {p},; includes the normalizations 
of t he nearby point and diffuse sources and the J factor, 
J*. log 10 ( Jj) and a ,■ are the mean and standard devia- 
tions of the distribution of log 10 (Jj), approximated to be 
Gaussian, and their values are given in Columns 5 and 
6, respectively, of Table I. 

The fit proceeds as follows. For given fixed values of 
mw and bf, we optimize — InL, with L given in Eq. 1. 
Confidence intervals or upper limits, taking into account 
uncertainties in the nuisance parameters, are then com- 
puted using the ‘‘profile likelihood” technique, which is 
a standard method for treating nuisance parameters in 
likelihood analyses (see, e.g., [32]). and consists of calcu- 
lating the profile likelihood — In L p ((a ann v)) for several 
fixed masses mw , where, for each (<7 ann u), — Ink is min- 
imized with respect to all other parameters. The inter- 
vals are then obtained by requiring 2Aln(£ p ) = 2,71 for 
a one-sided 95% confidence level. The MINI 11' subrou- 
tine MINOS [33] is used as the implementation of this 
technique. Note that uncertainties in the background fit 
(diffuse and nearby sources) are also treated in this way. 
To summarize, the free parameters of the fit are (a ana v ) , 
the J factors, and the Galactic diffuse and isotropic back- 
ground normalizations as well as the normalizations of 
near-by point sources. The coverage of this profile joint 
likelihood method for calculating confidence intervals has 
been verified using toy Monte Carlo calculations for a 
Poisson process with known background and Ferrm-LAT 
simulations of Galactic and isotropic diffuse gamma-ray 
emission. The parameter range for (cr tinn v) is restricted 
to have a lower bound of zero, to facilitate convergence of 
the MINOS fit, resulting in slight overcoverage for small 
signals, i.e., conservative limits. 


RESULTS AND CONCLUSIONS 

As no significant signal is found, we report upper lim- 
its. Individual and combined upper limits on the anni- 
hilation cross section for the bb final state are shown in 
Fig. 1: see also [34]. Including the J-faetor uncertainties 


Upper limits, bb channel 


... 3-10-* 

- - Draco 

- - Sextans 

- Bootes I 

Fornax 

- - Ursa Major II 

- - Carina 

— ~ Sculptor 

— - Ursa Minor 

— Coma Berenices 

- — Segue 1 

— Joint Likelihood, 10 dSphs 



io 1 to 1 to 3 

WIMP mass [GeV] 


FIG. 1. Derived 95% C.L. upper limits on a WIMP anni- 
hilation cross section for all selected dSphs and for the joint 
likelihood analysis for annihilation into the bb final state. The 
most generic cross section (~ 3 ■ 10 26 cm 3 s -1 for a purely s- 
wave cross section) is plotted as a reference. Uncertainties in 
the J factor are included. 



FIG. 2. Derived 95% C.L. upper limits on a WIMP annihila- 
tion cross section for the bb channel, the t + t~~ channel, the 
p + p _ channel, and the W~W~ channel. The most generic 
cross section (~ 3- 10~ 26 cm 3 s“ 1 for a purely s-wave cross sec- 
tion) is plotted as a reference. Uncertainties in the J factor 
are included. 


in the fit results in increased upper limits compared to 
using the nominal J factors. Averaged over the WIMP 
masses, the upper limits increase by a factor up to 12 
for Segue 1, and down to 1.2 for Draco. Combining the 
dSphs yields a much milder overall increase of the upper 
limit compared to using nominal J factors, a factor of 
1.3. 

The combined upper limit curve shown in Fig. 1 in- 
cludes Segue 1 and Ursa Major II, two ultrafaint satel- 
lites with small kinematic data sets and relatively large 
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uncertainties on their J factors. Conservatively, exclud- 
ing these objects from the analysis results in an increase 
in the upper limit by a factor ~1.5, which illustrates the 
robustness of the combined fit. 

We recalculated our combined limits using, for the clas- 
sical dwarfs, the J factors presented in [35], which allow 
for shallower profiles than Navarro-Frenk- White assumed 
here. The final constraint agrees with the limit from our 
J factors to about 10%, demonstrating the insensitivity of 
the combined limits to the assumed dark matter density 
profile. 

Finally, Fig. 2 shows the combined limits for all stud- 
ied channels. The WIMP masses range from 5 GeV to 
1 TeV, except for the W + W~ channel, where the lower 
bound is 100 GeV. For the first time, using gamma rays, 
we are able to rule out models -with the most generic cross 
section (~ 3 ■ 10 -26 cm 3 s -1 for a purely s-wave cross sec- 
tion), without assuming additional astrophysical or par- 
ticle physics boost factors. For large dark matter masses 
(around or above a TeV), the radiation of soft electro- 
weak bosons leads to additional gamma rays in the en- 
ergy range of relevance for the present analysis (see, e.g., 
[36, 37]). This emission mechanism is not included in the 
Monte Carlo simulations for the photon yield we employ 
here. While massive gauge boson radiation is virtually 
irrelevant for masses below 100 GeV, our results for the 
heaviest masses can be instead viewed as marginally more- 
conservative than with the inclusion of radiative electro- 
weak corrections. 

In conclusion, we have presented a new analysis of the 
Fermi-LAT data that for the first time combines mul- 
tiple (10) Milky Way satellite galaxies in a single joint 
likelihood fit and includes the effects of uncertainties in 
updated J factors, yielding a more robust upper limit 
curve in the (%,(% m v)) plane. This procedure allows 
us to rule out WIMP annihilation, with cross sections 
predicted by the most generic cosmological calculation 
up to a mass of ~ 27 GeV for the bb channel and up 
to a mass of ~ 37 GeV for the t + t~ channel. Future 
improvements planned by the Fermi-LAT Collaboration 
(apart from an increased amount of data) will include an 
improved event selection with a larger effective area and 
photon energy range and the inclusion of more satellite 
galaxies. 
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Note added in proof.- During the final preparation for 
submission of this Letter, we became aware of the work 
by Geringer-Sameth and Koushiappas when it was posted 
to the arXiv [38], reaching similar conclusions as ours, 
albeit with a different analysis. 
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